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Holstein-Hubbard model at half filling : A static auxiliary field study 


Saurabh PradharH and G. Venketeswara Pai 
Harish-Chandra Research Institute, Chhatnag Road, Jhunsi, Allahabad 211019, India 

We study the Holstein-Hubbard model at half hlling to explore the ordered phases such as the 
charge density wave and antiferromagnet. The Coulomb interaction is rewritten in terms of auxiliary 
helds. By treating the auxiliary fields and phonons as classical, we obtain real space features of the 
system and transition between the phases from weak to strong coupling. When both interactions 
are weak, mutual competition between them leads to a metallic phase in an otherwise insulator 
dominated phase diagram. Spatial correlations induced by thermal fluctuations lead to pseudogap 
features at intermediate range of coupling. 
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I. INTRODUCTION 

Several materials, notably transition metal oxideJ^^, 
have strong Coulomb interactions among their con¬ 
stituent electron^, as well as strong coupling between 
electrons and the underlying lattic^. Interplay of 
such competing many body interactions often leads to 
the emergence of effective energy scales and various 
broken symmetry phases and transitions among them, 
giving r ise to significant changes in their low-energy 
behavioil^. Understanding the combined effect of 
these two is a challenging problem and there has 
been some remarkable progress in the last couple of 
decades. Some of the systems that are known to 
have both these interactions playing a role include 
high-Tc superconducting cupra tes^t^, alkali doped 
fullerideJi^HlIl^ b ismut hate j^^ft^ j and most notably 
doped manganite d^^ l ^^ l In cuprated^*^, kinks observed 
in ARPES are believed to be features arising from 
strong electron-phonon coupling which also give rise to 
prominent features in inelastic neutron scattering and 
tunneling. The system also has strong electron-electron 
interactions as evidenced by the Mott insulating state 
of the parent compound. In fullerides, an antiferro¬ 
magnetic phase stabilized by Coulomb interactions, 
evolves to an s-wave superconducting state and it is 
believed that phonon effects are likely to be present. 
In doped bismuthates, a charge density wave^^ trans¬ 
forms to an s-wave superconducting state upon doping; 
valence skipping arising due to Coulomb interactions 
and coupling of charge carriers to breathing mode 
phonons are believed to be responsible for the behavior. 
Manganited^sHU] present the most compelling case 
where orbitally degenerate electrons experience strong 
Mott-Hubbard interactions and are also coupled to 
octahedral symmetry lifting Jahn-Teller phonon modes. 
It is being realized that the conventional way of treating 
only one of the interactions is inadequate for a proper 
understanding of these materials. 

The Holstein-Hubbard modePHm is the simplest 
starting point to theoretically explore the combined ef¬ 


fect of these two interactions. It describes a single-band 
electron coupled to an Einstein phonon mode. The 
Coulomb interaction is modeled by an on-site Hubbard 
term capturing the energy cost when two electrons of 
opposite spins are present at a given site. In real sys¬ 
tems, this model could be an oversimplification. There 
could be multiple orbitals relevanlP^ as in manganites, 
leading to inter- and intra-Coulomb matrix elements. 
There could also be multiple phonon modes involved 
as happens in Jahn-Teller system^^. However, general 
features, leaving out specifics such as orbital ordering, 
would be very well captured by the simplest model itself. 
Eor example, on a two-dimensional square lattice at 
half filling, the Hubbard interaction is expected to give 
rise to a weak-coupling spin density wave transforming 
to a local-moment antiferromagnetic Mott-Hubbard 
insulator (MHI) at strong coupling accompanied by a 
metal-insulator transition at finite temperature^^IIlll 
On the contrary, the Holstein interaction promotes 
coexisting charge density wave and superconducting 
ground states, if phonon dynamics is retained. However, 
for static phonons, it is expected that a weak coupling 
charge density wave would crossover to a bipolaronic 
insulator at strong coupling. Obviously, these phases 
will compete strongly when both interactions are 
present. Motivated by this, there have been several 
studies in recent years. These include analysis of 
various aspects of the problem using Migd al-Eli ashberg 
theorjl^, quantum Monte Carlo (QMCj^fiUSI^ exact 
diagonalisatioifil*^, variational treatment^^ such as 
the Gutzwiller approximatiorP^ for correlat ion, and 
dynamic mean field theory (DMFTj(23l2311HlSI, jn 
particular, Bauer and HewsoiP^ studied the grou nd 
state of the model at half filling using OMFl^^^^HIS] 
in co njunct ion with numerical re-normalization group 
(NRGp^*^. A recent studj® using dynamical mean 
field theory with continuous-time quantum Monte Carlo 
as an impurity solver has brought out several interesting 
features . These include strong renormalization of su¬ 
perconducting Tc and the emergence of a paramagnetic 
metallic phase in the weak-coupling limit. While DMFT 
is by far one of the most reliable tools to study strongly 
correlated systems, it has certain limitations. It is exact 
in infinite dimensions or when coordination number is 
large; however the theory being a local one does not 
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capture the full real-space features. If the system has 
geometrical constraints or frustration, a local theory 
will not be able to shed light on features intrinsic to 
them. It is also not possible to include disorder in any 
meaningful way since crucial features of interference 
cannot be captured in a single-site theory. Some new 
techniques are needed to overcome these problems 
and complement DMFT in the cases mentioned above. 
While such methods too will have their own limitations 
and range of applicability, they may be able to explore 
systems that DMFT cannot handle, especially when 
they are bench-marked in known cases. We use such a 
method to explore the problem at hand. 

The methocP^I^^IllHlIl includes rewriting the quar- 
tic fermion interaction in terms of auxiliary helds 
corresponding to charge and spin degrees of freedom. 
However, the resulting problem is still a many-body 
one, albeit with new fields. To simplify m atters , we 
concentrate on the static part of these held^SHISil 
also assume that the phonons are static. This results 
in a problem of a single-band electron moving in the 
background of three classical fields: the charge and 
magnetization auxiliary fields and lattice displacements. 
At any temperature, the statistically significant conhg- 
urations of classical fields ca n be sampled employing a 
Monte Carlo (MC) procedur d^ l ^^ l The electron problem 
can be solved by exact diagonalization. The method 
captures both weak and strong coupling regimes as 
described in Sec. III. 

We find that when only the Hubbard interaction is 
present, the system evolves from a Slatei!^ to a Mott 
insulator (MI^ with nonmonotonic variation of the Neel 
temperature. When only the lattice coupling is present, 
it transforms from a weak-coupling charge density 
wavdi^ to a bipolaronic insulator at strong coupling. 
When both are present, a critical line separates the two 
phases. At finite temperatures, the disordered phase 
appears to be metallic at weak coupling, but insulating 
at strong coupling. However, at intermediate coupling 
significant pseudogap features appeai^^ in the spectral 
function that modifies response of the electronic system 
such as optical transport in a significant way 

The paper is organized as follows. In Sec. H, we 
describe the model in detail and the method employed. 
Section HI is devoted to benchmarking with previous 
studies when only one of the interactions is present. In 
the next section, we give the ground-state (low temper¬ 
ature) phase diagram of the model resulting from the 
present study, followed by a detailed finite-temperature 
analysis of the electronic properties. Finally we conclude, 
describing the limitations of the method, advantages it 
has, and spell out future plans. 


II. MODEL AND THE STATIC AUXILIARY 
FIELD METHOD 

As mentioned earlier, we look at the simplest model 
of a one-band electronic model coupled to a single-mode 
Einstein phonon with the Coulomb interaction assumed 
to be local on a two-dimensional square lattice. The 
Hamiltonian is given by 

H T ^el—ph “b J^ph: 

Htb = ~t 'y ) cl^Cja- + H.C., 

(bl.o- 

^Hubbard — U ^ 
i 

H+ y H Q"’ 
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^el—ph — 9 ■ (^- 1 ) 
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Here Htb is the kinetic energy of the electronic system 
with t being the hopping parameter, Ci being an electron 
destruction operator at site f, and (ij) representing the 
nearest neighbors j of site i. U is the on-site Hubbard 
interaction and g is the Holstein electron-phonon cou¬ 
pling. Hpfi is the Hamiltonian for the Einstein phonon 
with frequency oo = ^K/m. Since we are interested 
in half filling {rii) = 1 . Since the classical single-site 
Holstein Hamiltonian has a polaronic minimum with 
a distortion p = {g/K), and polaronic binding energy 
Epoi = — (^g^l2K), we scale the phonon coordinate Q 
by p and phonon energies by \Epoi\. This results in 
a single dimensionless parameter (scaled in terms of 
energy unit t) for the phonon part of the Hamiltonian 
which we denote as V. From now on, we denote the 
dimensionless Hubbard interaction (in units of t) as U. 
We shall explore the physics of this model as functions 
of these two dimensionless parameters. 

To simplify this many-bo dy pr oblem, we perform a 
Hubbard-Stratanovich (HS)2SE21 transformation of the 
quartic interaction term by introducing two auxiliary 
fields, one each for the charge and magnetization 
sectors,. The scalar-valued charge auxiliary held at 
each site is (^i(T) and the vector-valued magnetization 
auxiliary held is mi(T). Since nt^nii = n?/4 — (s^ • mi) , 
where i J 2 a,i 3 cla^apCip, we can write 


g(7nit”i4- 


dfjjidmi 



-b idiTii + 


U 


- 2m, 



( 2 - 2 ) 

This results in a quadratic fermion problem in which 
fermions move around in a (quantum-mechanical, time- 
dependent) background of the two auxiliary helds and the 
phonon held which is computationally, still, a challenging 
problem. The partitions function is given by 
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+ i(j)^ni + - 2m* • 


(2.3) 


where Cph is the phonon Lagrangian. 


We make the following approximations. We assume 
that all three background fields are classical and hence 
neglect their time dependence. We retain their spatial 
dependence and do a thermal averaging of their configu¬ 
rations at every temperature numerically. We limit our¬ 
selves to half filling, i.e., one electron per site, in this 
paper. In this spirit, we make a saddle-point approx¬ 
imation for the static charge field, i.e., (j)i —> (0) = 
(C//2)(n*) = 17/2, and this is taken to be site inde¬ 
pendent. Upon res caling nii —)• ([7/2) nii, the resulting 
Hamiltonian reads jSSEII 

Heff = c|j,.Cjo- — ^effN — — nii • di 

(b>.o’ * 

+ ^Y.m^ + VJ2n^Q^ + vY,Ql (2.4) 

i i i 


The probability distribution functions appearing above 
are not exactly calculable since they involve tracing over 
fermions and integrating over all static configurations 
of the classical fields. We generate the equilibrium con¬ 
figurations for the classical field self-consistently using 
a Monte Carlo methocP^. This is achieved by starting 
with a given set of configurations, and attempting an 
update which requires diagonalizing the fermion Hamil¬ 
tonian and generating most probable configurations 
using the standard MC method. However, this severely 
restricts the system size of the problem, even though the 
fermionic part is quadratic. To explore higher system 
sizes, we use a traveling cluster algorithni^, in which 
a small cluster around the reference site is diagonlized 
and energy cost evaluated for MC update. During the 
MC procedure, as the reference site keeps moving on 
the lattice, so does the cluster. The results presented in 
this paper employ a cluster size of 8 x 8 and the largest 
system size used is 32 x 32. Once the system reaches 
equilibrium, we evaluate thermal averages of structure 
factor for charge density and magnetization. 

= (2-8) 

ij 

Spectral and transport properties for the fermion system 
have also been evaluated in thermal equilibrium which is 
described in Sec. V. 


where /ie// — /i — 17/2 with the partition function being j-jj_ EXPLORING THE HUBBARD AND 

given by HOLSTEIN PHYSICS 


Z = y DmX>QD[c^c] exp(-/3i7e//). (2.5) 


For a given configuration of Qi and mi, the Hamil¬ 
tonian (quadratic in fermions) needs to be diagonalized 
just once. However, one needs to sample most probable 
configurations of both Qi and at every temperature 
and they have to be determined from corresponding dis¬ 
tributions : 


P[Qi) 


P{mi) 


f VniD [ci, c] e f 
f Vm'DQ'D [ct, c] ’ 

JVQV [c\c] 
f 'Dm.VQ'D [ct, c] 


( 2 . 6 ) 

(2.7) 


While it appears that the neglect of the time- 
dependent effects reduces this method to unrestricted 
Hartree-Fock for the ground state, it retains the full 
classical thermal fluctuations in an unbiased way which 
leads to significant changes from HF results at finite 
temperature and smoothly interpolates between known 
limits at weak and strong coupling. 


In this section, we present the phase diagram of 
the model for the individual cases when either the 
Holstein term is absent (the Hubbard model) or the 
Hubbard term is absent (the Holstein model). Both 
these problems have been studied extensively in the past 
and it will help us benchmark our results. 

When the Holstein term is absent, the model reduces to a 
single-band Hubbard model on a two-dimensional square 
lattice at h alf fillin g. This does not have a metallic 
ground stat^^UllIllI for any nonzero value of U and has 
long-range antiferromagnetic insulating (AFI) order in 
the ground state. For small U, a Slater instability results 
due to nesting of the Fermi surface and the system is 
a spin density wave with a gap in the spectrum. For 
large U, the physics of superexchange takes over, due to 
the ”no double occupancy constraint” and the resulting 
kinetic energy reduction due to virtual hopping. The 
system is a Mott-Hubbard insulator with local moments 
present whose low-energy properties are governed by 
the antiferromagnetic Heisenberg model. The magnetic 
transitions resulting from these two behaviors have very 
different U dependence. At small U the Tn scales with 
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U as expected in an unrestricted HF treatment and 
results in a paramagnetic metallic phase (PM) above 
T^r due to the closing of the Slater gap . However, 
at large U, ^ i^/U) due to Mott physics and 
results in a paramagnetic insulator (PI). The present 
method captures both these behaviors very well. The 
finite-temperature phase diagram also looks qualitatively 
different from the HF phenomenology. While for small 
[/, the Slater gap closes at T/v, there is a pseudogap 
(PG) state that appears at intermediate values which 
crosses over to a paramagnetic Mott-Hubbard insulating 
state at large U. The paramagnetic state has strong AF 
fluctuations, especially in the intermediate range of the 
coupling constant, which results in pseudogap features 
in the spectral function. 

We now consider the case when the Hubbard term 
is absent, resulting in the half-filled Holstein model on a 
two-dimensional square lattice. This model again does 
not have a metallic ground state for any nonzero V. For 
small V, there is a Peierls instability leading to a charge 
density wave due to nesting, with a gap in the spectrum, 
which we call a charge-ordered insulator (COI). The 
charge modulation occurs at a wave vector (7r,7r). At fi¬ 
nite temperatures, the gap shrinks and vanishes at Tcdw 
above which the system is a nonmagnetic metal (NMM) 
with passive spin degrees of freedom. As V increases, 
the ground state of the system evolves through this 
charge-ordered state resulting in a bipolaronic insulator 
(BPI) at very large V. This can be understood due to a 
mechanism similar to superexchange for the spins. Since 
U is absent, there is no energy cost for double occupancy 
and a bipolaronic state lowers the energy through virtual 
fluctuations of charge. In this limit the physics can be 
described using a nearest-neighbor-interaction model, 
i.e., H — where a ^ (l/I^) and hence 

the charge-ordering temperature goes as {1/V). At 
intermediate values of V, a pseudogap phase intervenes 
which has spectral and transport features similar to 
the one previously mentioned. The finite-temperature, 
large-14 phase is insulating with charges remaining 
as bipolarons, but losing their long-range order.The 
spin degrees of freedom are passive in the entire phase 
diagram and the magnetization vanishes. We give the 
two phase diagrams in the U — T and V — T planes in 
Fig. 1. 

In passing, we wish to point out that the above 
results are indeed not exactly what is expected in two 
dimensions since there cannot be any finite-temperature 
transitions. These results should be taken as sugges¬ 
tive of what would happen in higher dimensions or 
as crossover scales where correlation lengths increase 
rapidly. (See the concluding section.) Further, we 
characterize the pseudogap phase as one in which the 
density of states does not have any perceptible hard 
gap, but has a dip at the chemical potential, suggesting 
a dramatic decrease of the low-energy spectral weight. 




FIG. 1. (Color online) (a) The phase diagram of the Hub¬ 
bard model on a two-dimensional square lattice at half fill¬ 
ing. AFI, PM, PI, PG represent antiferromagnetic insulator, 
paramagnetic metal, paramagnetic insulator and pseudogap 
phases, (b) The phase diagram of the Holstein model on a 
two-dimensional square lattice at half filling. COI, NMM rep¬ 
resent charge-ordered insulator and nonmagnetic metal. 


There is no real phase transition occurring here. It 
should be thought of as a crossover to a region where 
the density of state appears quite different from that of 
an insulator with a hard gap. 


IV. GROUND STATE PROPERTIES AND 
PHASE TRANSITIONS OF THE 
HOLSTEIN-HUBBARD MODEL 


Having clarified the trends that one obtains for the 
Hubbard and Holstein interactions separately, we now 
proceed to discuss the results for the full problem. 
However, in this section, we will concentrate on the 
ground-state properties and the nature of the phase 
transitions at finite temperatures. This includes the 
U — V phase diagram at T =0, spectral functions 
of the fermions, probability density functions for the 
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lattice variables, and charge and magnetization field 
configurations in real space. The above information 
would help us correlate various trends and elucidate 
the physics that emerges. As the phases change while 
changing parameters, we will see that correlated changes 
occur in properties of the fermionic, phononic, and 
auxiliary field variables. 

In Fig. (2) we present the ground-state phase dia¬ 
gram of the Holstein-Hubbard model as a function of 
U and V. As expected the results along the horizontal 
and vertical axes (corresponding to cases when one 
of the parameters is absent) confirm the discussion 
in the previous section. The phase diagram is almost 
entirely dominated by insulating regions. This is not 
surprising since individually, each interaction tries to 
localize electrons giving rise to a band /Mott/bipolaronic 
insulator. In the intermediat^^lllll to large values of 
the scaled parameters, there is a transition between a 
charge-disordered magnetic insulator to a charge-ordered 
nonmagnetic insulator. For example, at large values of 
U, an otherwise MI in the absence of Holstein interaction 
transforms to a BPI as V increases. This is a result 
of the two competing interactions. While a large U 
tries to localize individual electrons at every lattice 
site at half filling, the Holstein interaction develops 
bipolaronic instability as discussed in the last section. 
When the energy scales become comparable, the system 
develops an instability and moves from one to the 
other..Notice that the spin structure factor at (7r,7r) is 
nonzero in the AF phase and S'(q) = 0 in the BPI phase 
signaling a nonmagnetic state. Similarly, the charge 
structure factor has a peak at q = (0,0) in the AF phase, 
and the modulation vector changes to (tt, tt) in BPI 
phase. At intermediate values of U and V this behavior 
persists for both the structure factors but is much less 
pronounced compared to the strong-coupling limit. 
This is the crossover regime between the Slater-MHI 
due to Hubbard correlations and Peierls-BPI crossover 
due to Holstein interaction. Figure 3 depicts trends 
of both the structure factors as a function of U and 
V and confirms our conclusion about the transitions. 
Previous studied on the t — J-Holstein model have 
shown that as the exchange J decreases, the critical 
electron-phonon coupling required for the transition 
from AFI to COI increases. This is consistent with our 
results since J ~ \/U. However, there exists a thin 
sliver of window in the U — V plane at low interaction 
strengths where the system is metallic. This is in 
contrast to the case where the system is insulating when 
only one of the interactions is present. This behavior 
is exemplified in Fig. (3) where the structure factor at 
these values is plotted. This metallic behavior has been 
observ ed in previous studies of this model employing 
DMFT124121I 

using continuous-time QMC and NRG as 
impurity solvers. This unexpected m etallic phase results 
from the fact that while the nestin^i^*^ at half Filing 
in the two-dimensional tight-binding model supports 


magnetic or charge-ordering instabilities separately, the 
competing interactions have a destructive effect on the 
transition since it frustrates different degrees of freedom, 
viz., charge and spin in our case. The energy gained by 
a small mean field gap opening up in either channel is 
not sufficient to lower the absolute ground state energy 
when the other channel is included. This phase, in 
fact, brings out the true competition between the two 
interactions, where one acts predominantly over the 
spin sector while the other over the charge sector. A 
previous DMFT studj^H has revealed that this metallic 
region expands to larger U and V values as phonon 
frequency increases, which could explain the stability 
of this phase in the classical phonon limit. Further, 
notice that the phase boundaries merge to zero values 
of both parameters in our case in contrast to DMFT 
results. This is easily understood since our method 
preserves the nesting instability of the two-dimensional 
non-interacting electron system whereas methods such 
as DMFT ignore them. 


The MC procedure allows us to track the PDF of 
the phonon displacement variables across the transi¬ 
tions/crossover which is plotted in Fig. (4). In the AF 
phases we see that P(Q) is a unimodal function peaked 
at Q = -I, which implies that while every lattice site 
is distorted, it accommodates at most one electron per 
site. The distribution grows sharper as we grow from 
Slater to Mott limit, but the unimodal nature does 
not change. In this limit, Hubbard correlations play a 
larger role and the system tries to reduce the maximum 
number of electrons to one per site. At intermediate 
and strong coupling, at fixed U as we increase V, we 
find that this unimodal distribution slowly crosses over 
to a bimodal one. This occurs because of the weakening 
of the Hubbard correlation and increasing role of the 
polaronic distortion energies. Two electrons of opposite 
spin occupying the same site lower the electron phone 
energy more and the system develops a bipolaronic in¬ 
stability [also see Figs. (2a) and (2b)]. In the nontrivial 
metallic phase, while every site is still distorted, the 
amplitude is very small. These results, indeed, correlate 
with the charge structure factor and phonon probability 
distribution function. 

Our method allows us to provide a direct picture of 
the real-space correlations between the static magnetic 
auxiliary held, charge density, and phonon variables at 
various sites. This will elucidate the character of the 
transition and especially the metallic phase that arises. 
In Fig. (5) we present snapshots of spin and charge over 
the lattice for a given set of parameters at a given instant 
of MC simulation after the system has equilibrated. 
Nonlocal correlations among charge and magnetization 
helds can, in principle, be extracted from here. As 
expected, for lower values of V, spin correlations develop 
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FIG. 2. (Color online) The ground-state and finite- 
temperature phase diagram of the Holstein-Hubbard model as 
a function of scaled parameters U and V. AFI, COI, NMM 
represent antiferromagnetic insulating, charge-ordered insu¬ 
lating and nonmagnetic metal phases. The transition between 
AFI and COI is a weak first-order one. The temperatures are 
(a) T = 0.001 and (b) T = 0.050. 
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FIG. 3. (Color online) Structure factor corresponding to 
charge density wave A^(7r, 7r)(left) and the antiferromagnetic 
structure factor S(ti, tt)( right) versus temperature for U = 4.0 
and some representative values of V. 
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FIG. 4. (Color online) Phonon probability distribution for 
different values olU &tV = 2.0, T — 0.001. 


as U increases moving to a local moment value in the 
MI phase. Such spin correlations are absent in the COI 
phase. On the contrary, charge densities modulate as 
V is changed for a given U resulting in a bipolaronic 
state. In the corresponding phases spin modulation 
is negligible. In the metallic phase, both densities 
remain negligible on average, but there are fluctuations. 
The snapshots show a given configuration with some 
variation in the densities. However, an average over such 
configurations results in uniform charge density and 
negligible magnetization confirming that it is indeed a 
metallic phase. 


The various physical quantities that we have used 
to characterize the ground-state properties confirm 
the expected behavior and is reflected in such diverse 
variables as charge, magnetization, distribution of lattice 
displacements, and fermion spectral functions. The 
real-space picture gives a handle on how to correlate 
them. As will be discussed in the final section, this gives 
the added advantage of visualizing such changes in non¬ 
trivial geometries and especially on frustrated lattices, 
which is intractable or computationally expensive using 
other methods such as DMFT or its cluster variants. 


To conclude this section, within the static auxiliary 
field approximation of the decoupled HS fields that we 
have resorted to, we find a phase diagram that at low 
values of electron-phonon coupling crosses over from a 
Slater to MH insulator as U is increased, and a Peierls to 
BPI as V is increased for low values of Hubbard interac¬ 
tion. At intermediate to large values of coupling, there 
is a transition from antiferromagnetic MHI to a nonmag¬ 
netic charge-ordered or bipolaronic insulator. However, 
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FIG. 5. (Color online ) Charge rii (upper) and spin configu¬ 
ration Sj.So (lower) for U = 2.0. Temperature increases from 
left to right. Center column shows the conhguration near the 
Tc. Three rows for V — 0.50 (top), 1.0 (middle), 2.0 (bottom) 
and the system size is 32 x 32. 


there is a sliver of metallic phase at low coupling that re¬ 
sults from frustrating effects of two interactions in differ¬ 
ent (charge and spin) channels. The behavior of different 
degrees of freedom correlates with these changes provid¬ 
ing us with an efficient way to extract physics from weak 
to strong coupling. 


V. SPECTRAL AND TRANSPORT 
PROPERTIES 


Significant changes are expected in the phase diagram 
at finite temperatures due to the inclusion of ’’full” ther¬ 
mal fluctuations of the static field through configuration 
sampling. This was already noted in Sec. Ill where the 
effect of each interaction was looked at separately. In 
this section we present the results for various physical 
properties at finite temperatures and converge on the 



FIG. 6. (Color online) Density of state for different values 
temperature at constant V = 2.00 and U varying across the 
charge density wave-antiferromagnetic transition. U = 1.0 (a) 
, 4.0 (b), 6.0 (c), 8.0 (d). 


finite temperature phase diagram. 

Figure 6 shows the thermal-averaged single-electron 
spectral function A{uj) for different parameters at differ¬ 
ent temperatures. Deep in the insulating phase and at 
low temperatures they show a very clear gap and there 
are no states available at the Fermi energy as expected. 
In the region where metallic ground state appears, on the 
contrary, there is nonzero spectral weight at the Fermi 
energy even at the lowest temperatures. As temperature 
increases, we notice three regimes signifying different 
spectral features. For large values of U and/or V, we 
find that the gap persists even for large temperature. 
This is due to the Mott-Hubbard or bipolaronic nature 
of the phases. The fact that this feature survives at these 
values of parameters shows that the present method is 
capable of capturing the strong-coupling physics of this 
problem in both channels. At weak couplings, where 
a Slater or Peierls insulating phase is expected or the 
metallic phase emerges, the spectral features are very 
different at high temperatures. The gap vanishes entirely 
in the former cases and there is sufficient weight at the 
Fermi energy in all the three regimes. This clearly shows 
that the gap arises solely due to the nesting instability 
of the underlying Fermi system and the resulting order 
in either spin or charge channels. Once the order is 
destroyed, so is the gap. The most interesting features 
arise at intermediate values. Here a hard gap is not seen 
though there is significant reduction of spectral weight 
near the Fermi energy. There is spectral weight transfer 
from the coherence peaks to energies within the gap. 
This pseudogap feature arises due to persistence of local 
correlations in static fields even after the long-range 
order is destroyed. 


The MC snapshots throw more light on the existence 




















































































FIG. 7. (Color online) Temperatnre dependence of the optical 
conductivity for V — 2.0 and U = 1.0 (a), 4.0 (b), 6.0 (c), 8.0 
(d). 

of short-range order in either spin or charge degrees 
of freedom at temperatures near or above the ordering 
temperatures. This is shown in Fig. (5). In each case, 
the states evolve from the ground states shown in Fig. 
(5). However, unlike the low-coupling counterparts, the 
local order persists even above transition temperatures. 
This local order, we believe, is the reason for the appear¬ 
ance of pseudogap-like features in spectral functions. 
However, unlike the strong-coupling cases, where a local 
moment or a bipolaron formation is favored and the 
spectrum shows a hard gap, the intermediate range 
does allow fluctuations in charge and spin variables at 
very site, leading to spectral weight appearing in the 
otherwise gapped region. We have verified that the 
phonon PDFs also exhibit persistence of bimodality in 
this region. 

The transport can be captured in an exact way with¬ 
out resorting to approximations as in cluster DMFT. To 
this end, we use Kubo formulc!^ for the in-plane resis¬ 
tivity which involves the exact eigenvalues (ea,e/ 3 ) and 
wave functions (|a) , |/3)) of fermions obtained from diag- 
onalization at several equilibrium configurations. 

(") = V E u.i (-■ - 

^ - e. 

(5.1) 

Here / denotes the corresponding Fermi function and the 
current operator is given by 

Jx — it ^ cr h.c^ , (5.2) 

i,(T 

2 

where tJo = The dc conductivity is obtained by 

letting w —)• 0. 



T 


FIG. 8. (Golor online) Temperature dependence of resistivity, 
p(T), for various values of U at V — 2.0. Metal-insulator 
transition in the weak-coupling region for U = 1.0 (inset). 


Figure 7 shows the evolution of optical conductivity 
for a fixed value of V, but for varying U at different 
temperatures. A notable feature is the non-Drude be¬ 
havior of cr(a;). Further, the pronounced low-frequency 
hump in the optical conductivity for small frequencies 
evolves into an inter-band Hubbard peak as U increases. 
A similar feature has been observed as we vary V where 
the Hubbard peak gets replaced by the higher energy 
bipolaronic peak. The non-Drude behavior emanates 
from the pseudogap nature of the electronic spectral 
function that originates from strong local charge/spin 
fluctuations as discussed earlier. The dc resistivity is 
plotted in Fig. (8) for a fixed value of V, but varying val¬ 
ues of U at different temperatures. A metal-to-insulator 
transition is clearly visible for weak-coupling regime 
{U = 1) in the inset. 

Finally, we present the finite-temperature phase dia¬ 
gram of the model in Fig. (9). The phases include 
AF or CO insulating phases at low temperatures ex¬ 
cept for the sliver of metallic phase discussed earlier, 
metallic nonmagnetic phases at weak couplings, Mott- 
Hubbard and bipolaronic insulating phases at large cou¬ 
plings, and the pseudogap phase at intermediate cou¬ 
pling. The high-temperature behavior from metallic to 
insulating is a crossover. Note that we characterize the 
finite-temperature metallic phase by sign of the tem¬ 
perature variation of the resistivity, dp/dT. It remains 
open as to how these instabilities would be affected due 
to quantum dynamics of the auxiliary field or phonons. 
However, the remarkable qualitative agreement with pre¬ 
vious DMFT studies suggests that the quantum dynam¬ 
ics of these fields may not be relevant for the regime 
we have concentrated on. Further, it appears that the 
current method may be used for geometries and systems 
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FIG. 9. (Color online) Transition temperature for the charge 
density wave and antiferromagnetic phase for different values 
of U and V. 

where DMFT treatment may not be applicable as we dis¬ 
cuss below. 


VI. CONCLUSIONS 

We presented above a numerical study of the static 
Holstein-Hubbard model by employing Hubbard- 
Stratanovich auxiliary fields for the charge and spin 
sectors. It captures many of the features obtained in 
previous DMF studies. More importantly, it sheds light 
on new physics at finite temperatures and intermediate 
couplings due to the inclusion of spatial dependence 
(unlike DMFT) and classical thermal fluctuations 
through configuration sampling. The method works 
very well at all strengths of coupling. Being a real-space 
method it allows one to visualize the phases at various 
temperatures and how different orders develop and 
transform into others. It is numerically more efficient 
and large system sizes can be accessed. Various physical 
properties such as single particle spectral functions, 
phonon distributions, and transport can be readily 
evaluated. 

The salient results include the appearance of a 
nonmagnetic metallic phase at low values of coupling 
parameters, which was also seen in previous DMFT 
study, in addition to the ordered phases, including 
antiferromagnetic and charge-modulated ones. However, 
the finite-temperature phase diagram shows rich features 
and includes a pseudogap phase at intermediate cou¬ 
pling. This arises due to the persistence of local order in 
charge and spin degrees. Inclusion of spatial correlations 
is essential to capture this region. However, a couple of 
remarks are in order. First, the ground-state transition 
from AF to CO is expected to be a first-order one. How¬ 
ever, we find that this is a very weak transition and we 
are not able to resolve it accurately within the numerical 


error bars. The weak nature of this transition was also 
noted in previous DMFT studies. The low-temperature 
insulating states at small couplings could be a result of 
the fact that we used a two-dimensional square lattice. 
This necessarily gives a nesting instability at half 
filling and results in Slater or Peirels transition at low 
temperatures. Inclusion of quantum fluctuations or use 
of different lattice geometries may obscure these phases. 
Indeed, DMFT study shows that the metallic phase 
exists at low strengths even when one of them is zero. 
Third, one could wonder whether these transitions are 
numerical artifacts since we have used a two-dimensional 
system. Since thermal fluctuations destroy any order at 
nonzero temperatures, we expect Tjv = 0 and Tcdw = 0. 
However, it is expected that there would be a coherence 
temperature roughly mimicking the above transition 
temperatures even in two dimensions below which the 
correlation lengths increase rapidly. In other words, the 
system enters the renormalized classical regim^^. If so, 
even a weak coupling to a third dimension will stabilize 
the ordered phases. There could be some qualitative 
changes such as disappearance of insulating phases at 
weak couplings since nesting is no longer possible, but 
we expect gross features to remain the same. 

The method presented neglects time dependence in 
auxiliary fields and phonons. Comparison of our 
results with previous DMFT studies suggests that 
quantum dynamical effects may not be highly relevant 
for these phases especially since the system orders at 
low temperatures. However, this is indeed a handicap 
and does not allow us to explore other instabilities 
such as superconductivity. The present method may 
be thought of in the same spirit as spin wave theory 
applied to spin systems, whereby one starts with a 
classical ground state configuration for the spins and 
builds up quantum corrections perturbatively in some 
small parameter. The present method is a step towards 
implementing such a procedure for many constituents 
interacting among themselves. However, it goes beyond 
the analogy of spin wave theory in several aspects. We 
do not need to assume a classical ground state; the 
Monte Carlo procedure selects it naturally. The latter 
also incorporates thermal fluctuations. The procedure 
handles the weak to strong coupling limits in a unified 
way by smoothly interpolating between them. A further 
approximation is made by treating the charge fields 
at the saddle-point level. The charge fields {(p) turn 
out to be purely imaginary and their contribution to 
the classical action does not have a lower bound. This 
makes the classical Monte Carlo sampling of these 
fields very unstable and necessitates the saddle-point 
approximation which we have resorted to. The results 
obtained this way, for the pure Hubbard and Holstein 
models, respectively (see Fig. (1)), suggest that the 
saddle-point approximation does not affect the results 
severely. Finally, since quantum fluctuations of the 
auxiliary fields are neglected, the results become less 
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reliable at low temperatures and lower dimensions, 
where quantum effects may either destroy or modify the 
classical ground state as happens in spin systems. A 
natural way of capturing corrections would be to allow 
small amplitude fluctuation of the classical variables 
around their equilibrium value at a Gaussian level and 
look at the stability of phases. That necessitates a new 
line of study and we postpone it for the futur^^. 

The current work is not aimed at probing the physics of 
various families of experimentally accessible systems in 
which phonon dynamics is crucial along with electron 
correlation effects. Instead, we attempted to study a 
model system using some physically appealing simpli¬ 
fications. However, we believe that the results would 
provide us some insight into the behavior of systems in 
which relevant phonon frequencies are very large. Let 
us look at some rele vant experimental systems. For 
the cuprate familjEm, t ^ o.25 eV, U ^ 3 eV, a; ~ 75 
meV, and the dimensionless electron-phonon coupling 
constant H ~ 1. Thus the scaled Hubbard interaction 
is roughly three times the electron-phonon coupling 
and at half filling they would be Mott insulators with 
no charge-ordering tendencies. However, the phonon 
dynamics is relevant since the relevant frequencies are 
not very large. For fullerides, bandwidth is roughly 
0.6 eV, U ~ leV, w ~ 90 meV, and V lies in the 
range of 0.5-11^^^. Here again. Coulomb interactions 
dominate, though phonon dynamics is crucial. The 
family of ma nganites presents a more complex and rich 
scenaricP^I^. The parameter values relevant for these 
systems turn out to be t ~ 0.2-0.4 eV, U ~ 3-4 eV, 
uj ~ 0.05 eV (the energy of the Jahn-Teller phonons), 
while the Jahn-Teller polaronic energy Ejt 0.5-1 
eV. While the scaled electron-phonon coupling V is 
large (in the range of 1-4), the adiabaticity parameter 
7 = huj/Ep ~ 0.2-0.3, where Ep is the Fermi energy, 
and hence small. Naturally these are candidates more 
relevant for the present study, albeit at half filling. 
However, the physics gets more complicated due to a 
variety of reasons: multiorbital and multi-(JT)-phonon 


effects and the Hund’s rule coupling, not to mention 
the cooperative nature of the JT modes and associ¬ 
ated charge-ordering tendencies. A generalization of our 
method should be appropriate to studying these systems. 

The method can be expanded to study many prob¬ 
lems of current interest. We mention a couple of them. 
In the present study we have limited ourselves to a 
single-band Hubbard model coupled to a single-phonon 
mode. Many interesting realistic systems, such as 
manganite d^^bs i jj-jdate^^, involve multiorbitals 

and multiphonon modes. However, the present method 
can be generalized naturally to include them. The 
computational complexity increases marginally, but the 
problem is tractable within the approximations used. 
The method could also be extended to study interfaced 
and/or heterostructured^ of correlated, electron-phonon 
problems which are difficult to handle in conventional 
methods that are being currently used. 

A central feature of the method lies in capturing 
spatial correlations. This is essential for capturing 
features arising due to local order. More importantly, 
methods such as DMFT that neglect spatiul depen¬ 
dence are not suited to study geometried^^^^ where 
spatial features are important. A relevant case is the 
Holstein-Hubbard physics in frustrated geometries. We 
are currently pursuing this problem which shows rich 
physics including transition from charge-ordered phases 
to charge stripes, nontrivial spin orders, etc. These 
results will be presented elsewherd^. 
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